The close interaction between hypoxia-related proteins and metastasis in pancarcinomas

Many primary-tumor subregions exhibit low levels of molecular oxygen and restricted access to nutrients due to poor vascularization in the tissue, phenomenon known as hypoxia. Hypoxic tumors are able to regulate the expression of certain genes and signaling molecules in the microenvironment that shift it towards a more aggressive phenotype. The transcriptional landscape of the tumor favors malignant transformation of neighboring cells and their migration to distant sites. Herein, we focused on identifying key proteins that participate in the signaling crossroads between hypoxic environment and metastasis progression that remain poorly defined. To shed light on these mechanisms, we performed an integrated multi-omics analysis encompassing genomic/transcriptomic alterations of hypoxia-related genes and Buffa hypoxia scores across 17 pancarcinomas taken from the PanCancer Atlas project from The Cancer Genome Atlas consortium, protein–protein interactome network, shortest paths from hypoxia-related proteins to metastatic and angiogenic phenotypes, and drugs involved in current clinical trials to treat the metastatic disease. As results, we identified 30 hypoxia-related proteins highly involved in metastasis and angiogenesis. This set of proteins, validated with the MSK-MET Project, could represent key targets for developing therapies. The upregulation of mRNA was the most prevalent alteration in all cancer types. The highest frequencies of genomic/transcriptomic alterations and hypoxia score belonged to tumor stage 4 and positive metastatic status in all pancarcinomas. The most significantly associated signaling pathways were HIF-1, PI3K-Akt, thyroid hormone, ErbB, FoxO, mTOR, insulin, MAPK, Ras, AMPK, and VEGF. The interactome network revealed high-confidence interactions among hypoxic and metastatic proteins. The analysis of shortest paths revealed several ways to spread metastasis and angiogenesis from hypoxic proteins. Lastly, we identified 23 drugs enrolled in clinical trials focused on metastatic disease treatment. Six of them were involved in advanced-stage clinical trials: aflibercept, bevacizumab, cetuximab, erlotinib, ipatasertib, and panitumumab.

The genes with the highest mean frequency per genomic/transcriptomic alteration were identified as it follows: mRNA upregulation in COPS5 (0.  Table S1). Additionally, the alteration with the highest mean frequency was mRNA upregulation (0.047), followed by CNV amplification (0.010), putative driver mutation (0.0033), CNV deep deletion (0.0032), mRNA downregulation (0.003), and fusion gene (0.0004). We performed the Bonferroni correction for multiple comparison test to identify statistically significant alterations through the PCA-TCGA. Consequently, we found that mRNA upregulation and CNV amplification were statistically significant altered (P < 0.001) across different types of alterations ( Fig. 2A).  Tumor stages and metastatic status. Figure 3A shows the mean frequency of genomic/transcriptomic alterations across 17 pancarcinomas per tumor stage. The tumor stages with the highest mean frequencies of alterations were T4 (0.072) and T2 (0.072), followed by T3 (0.068) and T1 (0.060) (Table S19 to S22). However, the Bonferroni correction test for multiple comparisons did not show significant differences (P > 0.05) between tumor stages. Figure 3B shows the association between alterations and the tumor stage with its highest mean frequency. T2 stage presented the highest mean frequency of CNV amplifications (0.012), CNV deep deletions (0.004), and fusion genes (0.0004); T3 stage presented the highest mean frequency of mRNA downregulation (0.004); and T4 stage presented the highest mean frequency of mRNA upregulation (0.050), and putative driver   Figure 3C shows the mean frequency of alterations per metastatic status. The mean frequency of M0 status was 0.069 (n = 4254), meanwhile, M1 status was 0.072 (n = 263) (Table S23). Finally, the Mann-Whitney U test showed that differences in alterations between M0 and M1 were statistically significant (P < 0.001).
Hypoxia score. Buffa et al. used information from in vitro experiments combined with in vivo co-expression patterns regarding hypoxia-regulated genes and pathways in order to construct the Buffa hypoxia score 58 . That is, higher mRNA abundance signatures indicate higher levels of hypoxia. Therefore, we quantified tumor hypoxia levels of 5,249 individuals with 13 pancarcinomas using the mentioned score 58 . Figure 4A shows the Buffa HS mean per cancer type. The most hypoxic pancarcinoma was HNSC (29.6), followed by LUSC (26.8),  Protein-protein interactome network. The PPi network was performed to better understand the connectivity between hypoxia-related proteins and metastatic driver proteins using the String database and the Cytoscape software v.3.7.1 59,60 . Metastatic driver proteins were taken from the Human Cancer Metastasis Database (HCMDB), an integrated database designed to analyze large scale expression data of cancer metastasis 61 . Subsequently, we generated the interactome network encompassing 108 nodes and 603 high-confidence interactions (cutoff = 0.9) (Fig. 5). The mean of degree centrality of the PPi network was 11.2, and the top ten hypoxic proteins involved in metastasis with the highest degree centrality were AKT1 (43) Table S26. On the other hand, angiogenesis, the recruitment of new blood vessels, is an essential component of the metastatic pathway 63 . In our study we found 106 (45%) hypoxic-related proteins with shortest paths to angiogenesis, of which, 73 had positive regulation with an average distance of 3.5 and an average path length of 4.4, 20 proteins had negative regulation with an average distance of 3.9 and an average path length of 4.4, and 13 proteins with unknown regulation status had an average distance of 3.7 and an average path length of 4.7. The top ten hypoxia-related proteins with the shortest distance to angiogenesis were TGFA (0.9), TGFB1 (0.9), TGFB2 (0.9), TGFB3 (0.9), THBS1 (0.9), TIMP1 (0.9), TNF (0.9), VEGFC (0.9), VEGFA (1.2), and MMP2 (1.7) (Table S26). Figure 7A shows a Venn diagram encompassing the 73 most altered proteins (mean frequency > 0.068) from PCA-TCGA, 108 hypoxic/metastatic proteins with high-confidence interactions (cutoff = 0.9) from the PPi network, and 112 hypoxia-related proteins with the shortest paths to    64 . Figure 7B shows a circos plot encompassing 12 of the 30 essential hypoxia-related proteins already characterized as hallmarks of cancer according to the Catalogue of Somatic Mutations in Cancer (COSMIC) Cancer Gene Census (CGC) (https:// cancer. sanger. ac. uk/ census) 65 . As results, invasion and metastasis were promoted by EGFR, AKT1, MTOR, NOTCH1, ERBB2, MAPK1, and CDH1; escape from programmed cell death was encouraged by AKT1, EGFR, MTOR, NOTCH1, ERBB2, PIK3CA, and MAPK1; changes in cellular energetics were stimulated by TP53, AKT1, EGFR, MTOR, NOTCH1, ERBB2, ARNT, MAPK1, and CDH1; angiogenesis was endorsed by AKT1, EGFR, MTOR, NOTCH1, PIK3CA, and ARNT; proliferative signaling was supported by AKT1, EGFR, MTOR, NOTCH1, ERBB2, and PIK3CA; suppression of growth was promoted by TP53, AKT1, and EP300; escape from immune response during cancer was encouraged by EGFR; lastly, cell replication towards immortality was stimulated by TP53 and NOTCH1.

Functional enrichment analysis.
Subsequently, we performed a functional enrichment analysis of the 30 essential hypoxic/metastatic proteins obtained from the intersection between the PanCancer Atlas, the PPi network, and the shortest paths to metastatic and angiogenic phenotypes analyses. The functional enrichment analysis was performed through the g:Profiler software using the human cancer metastasis proteins as background set. We were able to identify 358 significant gene ontology (GO) biological processes, 99 significant the Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathways, and 65 significant Reactome signaling pathways [66][67][68] . The most significant GO biological processes with Benjamini-Hochberg correction and false discovery rate (FDR) < 0.001 were cellular response to hypoxia and positive regulation of cell migration. The most significant KEGG signaling pathways with Benjamini-Hochberg correction and FDR < 0.001 were HIF-1, PI3K-Akt, thyroid hormone, ErbB, FoxO, mTOR, insulin, MAPK, Ras, AMPK, and VEGF. Lastly, the most significant Reactome signaling pathways with Benjamini-Hochberg correction and FDR < 0.001 were PIP3 activates AKT signaling, PI3K/AKT signaling in cancer, signaling by ERBB2, and MTOR signaling ( Fig. 7C and Table S29). Drugs involved in clinical trials focused on metastatic disease. Figure 9 shows the current status of drugs in clinical trials for treatment of metastatic disease according to the Open Targets Platform 70 . There are 23 drugs targeting 10 hypoxic/metastatic proteins that have been analyzed in 211 clinical trials. The cancer types with clinical trials focused on metastasis were prostate cancer, colorectal cancer, and melanoma (Fig. 9A). EGFR was the target with the highest percentage of clinical trials in process, recruiting or completed, followed by ERBB2, VEGFA, IGFR1, NOTCH1, AKT1, AKT2, AKT3, ANGPT2, and PIK3CA (Fig. 9B). Most clinical trials were in phase 2 (62%), followed by phase 3 (19%), phase 1 (18%), and phase 4 (1%) (Fig. 9C). Tyrosine kinase EGFR family was the target class with the highest percentage of clinical trials (54%), followed by secreted protein (41%), tyrosine kinase InsR family (2%), AGC kinase AKT family (1%), and enzyme (1%) (Fig. 9D). Antibodies (77%) were the type of drugs focused on metastasis with the highest percentage of clinical trials, followed by proteins (13%) and small molecules (13%) (Fig. 9E). Lastly, Fig. 9F diagrams a Sankey plot comparing the number of clinical trials testing metastatic drugs in different cancer types. Bevacizumab (73), cetuximab (59), panitumumab (32), and aflibercept (14) were the metastatic drugs with the highest number of clinical trials (Table S30).

Discussion
Hypoxia seems to be a potential micro-environmental factor that induces metastasis through the activation of the HIF signaling pathway. The potential mechanisms evidenced include the reprogramming of metabolism, stem cell phenotype, invasion, vascular facts, suppression of immune system, pre-metastatic niche, intravasation, extravasation, and anti-apoptotic activity 71 . Poor vasculature induces imbalance between pro and anti-angiogenic signals, reducing oxygen delivery into tissues. The poor oxygen delivery produces low oxygen tension in tumor and stromal tissues 72 . In the process of tumor metastasis, invasion is observed to be the first step that is activated by pro-migratory factors induced by mesenchymal stem cells, collagen network formation, and recruitment of macrophages and fibroblasts 64,71 . Hypoxic tumor microenvironment is suggested to aggravate the metastatic initiation through the HIF signaling factor. Hypoxia also triggers the pro-migratory factors and activates the Scientific Reports | (2022) 12:11100 | https://doi.org/10.1038/s41598-022-15246-y www.nature.com/scientificreports/ extracellular matrix to support the metastatic process 42 . It is evidenced that the abnormal cancer cell intra and extravasation from the vascular structures is activated by the hypoxic tumor environment and the HIF factor which regulates the cell to cell endothelial adhesion molecules and promotes tumor metastasis 73 . Hypoxia induces a series of biological changes which contribute to tumorigenesis and the metastatic phenotype 28,74 , both of them are associated with resistance to radio-and chemotherapy, anti-cancer drugs, and immunotherapy. Thus, understanding the molecular signatures of hypoxia is crucial to identify potential therapeutic targets to improve metastatic disease treatments 75 . In this study, we revealed essential hypoxia-related proteins highly involved in metastatic signaling through three approaches that analyze genomic and transcriptomic alterations, protein-protein interactions, and shortest paths from hypoxia to metastatic and angiogenic phenotypes.
The first approach consisted in the analysis of 100,643 genomic/transcriptomic alterations (233 HRGs) belonging to 6343 individuals with 17 pancarcinomas. ESCA, BRCA and LUSC were the most altered cancer types, whereas MESO, KIRC, and THCA were the least altered. The mRNA upregulation was the most significant (Bonferroni correction, P < 0.001) genomic alteration since hypoxic tumors overexpress intracellular signals to adapt to the environment 16,76,77 . Overall, the first approach revealed that 73 HRGs presented frequencies of alteration higher than the average (> 0.068). It is important to mention that the 233 hypoxia-related genes were collected from gene ontology terms, signaling pathways, the Buffa signature, and publications on high-altitude adaptive phenotypes. Nevertheless, a limitation of this study is that we did not perform a manual search and curation of publications on individual genes potentially related to hypoxia in order to enrich our initial gene set.
Regarding to genomic/transcriptomic alterations per tumor stage across 17 pancarcinomas, T4 stage presented the highest mean frequencies of mRNA upregulation (0.050) and putative driver mutations (0.004). On the other hand, the mean frequency of genomic alterations in HRGs was significant in patients with metastasis (0.072) in comparison to patients without it (0.069) (Mann-Whitney U test, P < 0.001), meaning that the number of genomic alterations in the HRGs increases as tumor stage and metastatic status evolves. Smith et al.have proposed this statement through a genetic model in which sequential accumulation of mutations in specific genes (i.e., ACP, KRAS, SMAD2, SMAD4, and TP53) drives the transition from healthy epithelia to metastatic colorectal cancer 78 .
In order to validate the significant differences in alterations of HRGs found in tumor stages and metastatic status, we analyzed the Buffa hypoxia score across 13 pancarcinomas 58,79 . HNSC, LUSC, and CESC were the most hypoxic cancer types, whereas LIHC, PRAD, and THCA were the least hypoxic, as previously shown by Bhandari et al 16 . Regarding to tumor stages, T4 presented the highest HS mean (14.5) with significant differences (Bonferroni correction, P < 0.001) between T1 vs T4, T2 vs T4, and T3 vs T4. Regarding to metastatic status, the Buffa HS mean in patients with metastasis (8.3) was significantly higher (Bonferroni correction, P < 0.001) than in patients without it (2.2).  [80][81][82] , protein interactions with therapeutic significance can be revealed by the integration of cancer proteins into networks. Protein interactions regulate essential signals such as proliferation or metastasis, and thus, represent potential targets for drug development and drug discovery. Regarding our networking analysis, we generated an interactome network encompassing 108 nodes and 603 high-confidence edges (cutoff = 0.9) with a mean degree centrality of 11.2. Consequently, the second approach revealed 108 highly connected hypoxic and metastatic proteins. www.nature.com/scientificreports/ According to Ianuccelli et al., the value of the third approach permits to bridge the gap between genomic data and cancer phenotypes. However, it should be kept in mind that biology is more complex than graph theory. One important limitation of CancerGenNet is the finding that for up to 20% of the cancer genes we know little about the molecular mechanisms underlying their tumorigenic function 62 . In this context, the third approach was aimed at revealing the shortest paths from HRPs to metastasis 83 . Therefore, we found that 99 hypoxic proteins had paths leading to metastasis from which 49 had positive regulation and 16 had negative regulation. Additionally, the third approach was complemented with the analysis of angiogenic paths because the recruitment of new blood vessels is an essential component of the metastatic pathway 63 . Therefore, we identified 106 hypoxic proteins with shortest paths to the angiogenesis hallmark, of which, 73 had positive regulation and 20 had negative regulation. Then, the compendium of the 112 HRPs with the shortest paths to metastatic and angiogenic phenotypes made up the third approach.
Regarding clinical trials reported on the hypoxic/metastatic essential proteins, the Open Targets Platform is an available resource for the integration of proteogenomics and chemical data to aid systematic drug target prioritization 70 . In this study we identified 23 drugs targeting 10 hypoxic/metastatic proteins that have been analyzed in 211 clinical trials. Six of them (aflibercept, bevacizumab, cetuximab, erlotinib, ipatasertib, and panitumumab) were involved in phases III/IV clinical trials. Ipatasertib, erlotinib, lapatinib, neratinib, afatinib, sapitinib, gefitinib, tucatinib, linsitinib, and buparlisib are small molecules; aflibercept is a recombinant protein; and cetuximab, bevacizumab, panitumumab, futuximab, necitumumab, trastuzumab emtansine, pertuzumab, ganitumab, dalotuzumab, brontictuzumab, and MEDI-3617 are monoclonal antibodies. Lastly, the mechanism of action of these drugs is fully detailed in Table S31.
From a public health perspective, it is fundamental to recognize the holistic view that multi-omics provide to the understanding of cancer; for example, providing accurate selection of patients for the assessment of specific therapies. As reported by Atun et al., the cancer burden continues to grow globally, projecting to have around 24.6 million cases by 2030 85 . Physical health is not the only affected aspect, but millions of people who experience the economic, emotional, and social impacts of the disease and the increasing economic cost to the health systems. Specifically for low-and-middle income countries, there is a gap in the access to health, including a timely quality diagnosis and treatment. Additionally, several studies have highlighted the need to diversify oncological studies to populations representing several ethnic groups along with the development of novel strategies to enhance race/ethnicity data recording and reporting [86][87][88][89][90] . In this sense, multi-omics technologies and public-private investment related to identifying therapeutic targets improving metastatic disease treatments are crucial to reduce inequalities in health and strengthen mechanisms that can improve survival rates of different types of cancer. Relevantly, the effective and ethical use of these technologies would contribute to the development of knowledge and further explanations around the cause of diseases, with the ultimate aim of reducing morbidity and mortality; thus, increasing wellbeing. National health systems should work towards the identification of opportunity costs related to multi-omics research investment and its potential benefits on clinical applications to help diagnose and/or prevent certain diseases. As also stated by Hassin et al., multi-omics approach offers the opportunity to understand the flow of information that underlies disease 91 . The success of the omics approach has to be addressed not only in a technical and/or financial factors-assessment but on the importance of developing conceptual research shifts; focused also on the relations between biological factors and the interrelations with the social determinants of health.
In conclusion, hypoxia and HIF-dependent signaling play an important role in metastasis tumor progression. The hypoxic tumor microenvironment influences both the early and late stages of metastasis. Our findings suggest that individuals with metastasis present higher frequencies of genomic/transcriptomic alterations and Buffa hypoxia score across pancarcinomas from PCA-TCGA. The most altered signaling pathways were HIF-1, PI3K-Akt, thyroid hormone, ErbB, FoxO, mTOR, insulin, MAPK, Ras, AMPK, and VEGF. Finally, since cancer is a group of complex and heterogeneous diseases, the study of multi-omics approaches is an effective way to understand the molecular landscape behind hypoxia and metastasis that revealed 30 potential therapeutic targets and 23 drugs to improve metastatic disease treatments. These drugs can be considered for treating metastasis after being thoroughly evaluated in specific clinical trials.   55,56 . Lastly, the Mann-Whitney U test and the Bonferroni correction test for multiple comparisons were performed to determine significant differences (P < 0.001) between genomic/transcriptomic alterations, clinical data, and PCA-TCGA type.

OncoPrint of genomic and transcriptomic alterations.
Hypoxia score. The Buffa hypoxia score was analyzed in 5249 tumors from 13 pancarcinomas: HNSC, LUSC, CESC, CRC, BLCA, SKCM, KIRC, LUAD, PAAD, BRCA, LIHC, PRAD, and THCA 58 . An approach for deriving signatures that combine knowledge of gene function and analysis of in vivo co-expression patterns was used to define a common hypoxia signature in cancer. Hypoxia scores were estimated by obtaining the mean expression (log2) of 52 hypoxic genes reported by Buffa et al 58,79 . Data related to Buffa hypoxia score, tumor stage, and metastatic status were taken from the Genomics Data Commons and the cBioPortal 55,56 . Lastly, the Bonferroni correction test for multiple comparisons was performed to determine significant differences (P < 0.001) between hypoxia score and clinical data.
Protein-protein interactome network. The PPi network with high-confidence interactions (cutoff = 0.9) and zero node addition was created using the human proteome of Cytoscape StringApp 59,94 , which imports protein-protein interaction data from the String database 59,95-97 . The PPi network was encompassed by hypoxia-related proteins and metastatic driver proteins, which were taken from the Human Cancer Metastasis Database (https:// hcmdb.i-sanger. com/). HCMDB is an integrated database designed to analyze large scale expression data of cancer metastasis 61 . The degree centrality of a node represents the number of edges the node has to other nodes in the network and it was calculated using the CytoNCA app 98,99 . All nodes and edges were organized through the organic layout and visualized through the Cytoscape software v.3.7.1 60 . Lastly, the hypoxia-related proteins involved in metastasis signaling were differentiated by colors according to the degree centrality.
Shortest paths from hypoxia-related proteins to metastatic and angiogenic phenotypes. Can-cerGeneNet (https:// signor. uniro ma2. it/ Cance rGene Net/) is a resource that links proteins that are frequently mutated in all cancer types to cancer phenotypes. This resource is based on the annotation of experimental information that allows to embed the cancer proteins into the cell network of causal protein relationships curated in SIGNOR 100 . Therefore, this bioinformatics tool allows to infer likely paths of causal interactions linking cancer associated proteins to cancer phenotypes such as metastasis and angiogenesis 62,83,101 . Iannuccelli et al.explained that the shortest paths from a specific protein to cancer phenotypes was programmatically implemented using the shortest path function of igraph R package, obtaining a distance score and a path length score 83,102 . Hence, we analyzed the shortest paths from our hypoxic proteins to metastatic and angiogenic phenotypes to better understand the association to these hallmarks of cancer 62,64 .
Functional enrichment analysis. The enrichment analysis gives scientists curated interpretation of protein sets from omics-scale experiments 57,66,103 . Essential hypoxic and metastatic proteins were analyzed by using g:Profiler version e101_eg48_p14_baf17f0 (https:// biit. cs. ut. ee/ gprofi ler/ gost) to obtain significant annotations (Benjamini-Hochberg FDR < 0.001) related to GO biological processes, KEGG signaling pathways, and Reactome signaling pathways 66,68,92 . The functional enrichment analysis was performed using the Human Cancer Metastasis proteins as background set, and annotations were visualized through Manhattan plots. Lastly, significant terms related to hypoxia, metastasis, and oncogenic signaling pathways were manually curated.
Drugs involved in clinical trials focused on metastatic disease. The Open Targets Platform (https:// www. targe tvali dation. org) is comprehensive and robust data integration for access to and visualization of potential drug targets associated with several cancer types and metastasis. Additionally, this platform shows all drugs in clinical trials associated with hypoxic/metastatic proteins, detailing its phase, type of drug, action type, and target class 70 .

Data availability
All data generated or analyzed during this study are included in this published article (and its Supplementary Information files).